This research aims to explore different approaches (point-based analysis on reported damages and SAR imagery change detection) in performing a damage assessment that can guide the emergency response following the Beirut Port explosion that the capital on the 4th of August 2020.
Data included:
# loading libraries
library(spatstat)
library(here)
library(sp)
library(rgeos)
library(maptools)
library(GISTools)
library(tmap)
library(sf)
library(geojson)
library(geojsonio)
library(tmaptools)
library(dplyr)
library(stringr)
library(readr)
library(rgdal)
library(tmap)
library(janitor)
library(ggplot2)
library(raster)
library(fpc)
library(dbscan)
library(tidyverse)
library(tidyr)
# first, get the Beirut Explosion shapfiles
# operational zones reflects the damaged zones in Beirut
Bei_Exp_Zones <-
st_read(here::here("data", "beirut_port_explosions_operational_zones_139", "beirut_port_explosions_operational_zones_139.shp")) %>%
# transform the coordinate reference system to Dei EL Zor that covers Lebanon
st_transform(., 22770)
## Reading layer `beirut_port_explosions_operational_zones_139' from data source `C:\Users\SaraMoatti\Desktop\UCL\Year 2\GIS\Assignment\GIS_Final_Coursework\GIS_Final_Coursework\data\beirut_port_explosions_operational_zones_139\beirut_port_explosions_operational_zones_139.shp' using driver `ESRI Shapefile'
## Simple feature collection with 139 features and 15 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 35.49469 ymin: 33.86149 xmax: 35.5527 ymax: 33.90946
## geographic CRS: WGS 84
# get the beirut buildings geojson to reflect on the urban area
Beirut_Buildings <-
st_read("https://opendata.arcgis.com/datasets/d4d43fbe781145d4b11f9eac3f5dc5a1_0.geojson") %>%
st_transform(., 22770)
## Reading layer `Beirut_Buildings' from data source `https://opendata.arcgis.com/datasets/d4d43fbe781145d4b11f9eac3f5dc5a1_0.geojson' using driver `GeoJSON'
## Simple feature collection with 18340 features and 6 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 35.46775 ymin: 33.86234 xmax: 35.54173 ymax: 33.90645
## geographic CRS: WGS 84
# crop the dataframe with buildings to the operational zones extent
Bei_Zones_Buildings<- Beirut_Buildings[Bei_Exp_Zones,]
# get the UNHABITAT soscio-economic classification of operational zones
UN_habitat_shp <-
st_read(here::here("data","unhabitat_zone_data_batch2_2020_08_19", "UNHABITAT_Zone_Data_batch2_2020_08_19.shp")) %>%
st_transform(., 22770)
## Reading layer `UNHABITAT_Zone_Data_batch2_2020_08_19' from data source `C:\Users\SaraMoatti\Desktop\UCL\Year 2\GIS\Assignment\GIS_Final_Coursework\GIS_Final_Coursework\data\unhabitat_zone_data_batch2_2020_08_19\UNHABITAT_Zone_Data_batch2_2020_08_19.shp' using driver `ESRI Shapefile'
## Simple feature collection with 139 features and 14 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 35.49469 ymin: 33.86149 xmax: 35.5527 ymax: 33.90946
## geographic CRS: WGS 84
# rename the variable
UN_habitat_shp<-UN_habitat_shp %>%
dplyr::rename( socio_economic_classification = UNHABITAT1 )
#try to bring the rgb stack
#r1_stack<-
# stack(here::here("MAXAR_data","Raster tif before and #after the blast","Before","10300500A5F95600.tif"))
# unify the extent
# I won't do thi here as it takes ages due to the large raster
#r1_stack<-r1_stack %>%
# resample(.,Bei_Exp_Zones,
# method = "ngb")
# crop and mask
#r1_stack<- r1_stack%>%
# crop(.,Bei_Exp_Zones) %>%
# mask(.,Bei_Exp_Zones)
# repeat the process with image after the explosion
# after
#r1_after_stack<-
# stack(here::here("MAXAR_data","Raster tif before and #after the blast","After","10300500A6F8AA00.tif"))
# unify the extent
# I won't do thi here as it takes ages due to the large raster
#r1_after_stack<-r1_after_stack %>%
# resample(.,Bei_Exp_Zones,
# method = "ngb")
# crop and mask to Beirut operational zones
#r1_after_stack<- r1_after_stack%>%
# crop(.,Bei_Exp_Zones) %>%
# mask(.,Bei_Exp_Zones)
#check if both have same extent
#extent(r1_after_stack)
#extent(r1_stack)
#if not
#r1_after_stack<-r1_after_stack %>%
# resample(.,r1_stack,
# method = "ngb")
# start the image differencing with overlay
# trial with I2-I1
#stack_after_before <- overlay(r1_after_stack,
# r1_stack,
# fun = function(r1, r2) { return( r1 - r2) })
# save the outcome as it is easier to work with the image difference tif instead of repeating the lengthy process each time
#savetofolder_stack_after_before<-
# writeRaster(stack_after_before,
# filename=file.path(here::here("MAXAR_data","image_differencing"),"stack_after_before.tif"))
# increasing the memory to handle the process
#memory.limit()
#memory.limit(size=40000)
# log rationing difference
# log10(r1/r2)
#stack_after_before_log <- overlay(r1_after_stack,
# r1_stack,
# fun = function(r1, r2) { #return(log10(r1/r2))})
# save the outcome
#savetofolder_stack_after_before_log<-
# writeRaster(stack_after_before_log,
# filename=file.path(here::here("MAXAR_data","image_differencing"),"stack_after_before_log.tif"))
The code above won’t be run as files are very large and takes ages to process, this is to indicate how the data was prepared.
I will be calling the saved outcome : stack_after_before_log.tif to perform further analysis and save processing time but below is the outcome! let’s have a look at it!
stack_after_before_log<-
stack(here::here("MAXAR_data","image_differencing","stack_after_before_log.tif"))
# look at the result
plotRGB(stack_after_before_log,axes=TRUE, stretch="lin")
As the temporal difference is small, changes between before and after SAR images are assumed to be caused by the explosion, therefore, representing the damage.
Density Map per zone is generated, the raster::extract did not work as the file is huge the shapefile Raster_cell_count_per_zones.shp was created from QGIS through the zonal statistics plugin the image differencing tiff layer: stack_after_before_log and the beirut operational zones shapefile were used
# get the shapefile
Raster_cell_count_per_zones <-
st_read(here::here("MAXAR_data","image_differencing", "Raster_cell_count_per_zones.shp")) %>%
st_transform(., 22770)
## Reading layer `Raster_cell_count_per_zones' from data source `C:\Users\SaraMoatti\Desktop\UCL\Year 2\GIS\Assignment\GIS_Final_Coursework\GIS_Final_Coursework\MAXAR_data\image_differencing\Raster_cell_count_per_zones.shp' using driver `ESRI Shapefile'
## Simple feature collection with 139 features and 12 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 35.49469 ymin: 33.86149 xmax: 35.5527 ymax: 33.90946
## geographic CRS: WGS 84
# prepare the data
# get the density
density<-Raster_cell_count_per_zones%>%
#calculate area
mutate(area=st_area(.))%>%
#then density of the points per ward
mutate(density_raster=X_countfina/area)%>%
#select density and some other variables
dplyr::select(density_raster, zone_numbe, Cadaster_1, X_countfina)
tmap_mode("plot")
breaks = c(4, 4.15, 4.3, 4.45, 4.6,4.75)
dm1 <- tm_shape(density) +
tm_borders("white")+
tm_polygons("density_raster",
breaks=breaks,
palette="-cividis",
alpha = 0.7)+
tm_scale_bar(position=c(0.01,0.1),text.size=0.5,color.dark = "grey46")+
tm_compass(north=0, color.dark = "grey46",position=c("left","0.2"), size = 0.5)+
tm_legend(show=FALSE)+
tm_layout(title = "Damage assessment_SAR imagery based",title.size = 2,frame=FALSE)
dm1_legend <- tm_shape(density) +
tm_polygons("density_raster",
title = "Damage Density",
breaks=breaks,
palette="-cividis",
alpha = 0.7,
border.col = "white") +
tm_layout(legend.only=TRUE, legend.title.size = 1,legend.position=c(0.01,0.25),asp=0.1)
density_map1=tmap_arrange(dm1,dm1_legend)
density_map1
The below is an attempt to classify the raster cells/pixels where we can possibly understand the changes/damages
## this analysis follows the below tutorials
#http://remote-sensing.org/unsupervised-classification-with-r/
#https://www.r-exercises.com/2018/02/28/advanced-techniques-with-raster-data-part-1-unsupervised-classification/?__cf_chl_captcha_tk__=95093339d521789be882d0b5b3ab77f12debcf39-1609698000-0-AWqv_DC4nb9vdENGfsl9H-zfUjYkvB_MYG1F7eGd5_wbbLMYmiRCS7O4I-SO2gJ-vh0VHQjYruVihL_i6Yva7OXsHZSx6uy1ao-BCuTjPCDNdNI0rrL8rJZjC_-w-uTovyTYlZtZBqGyLMuMJAc1jkeYRoKKT-smd9quKJkufCXUuaKPOY0gOgKr3HEvnJ85RCLHlu5yR0nwwxE4Yby2OXixVxBC-zfs6KcIjSX_D0YPtFFcQJ6zeHLGvrKubUdgvy-QgK-yibPsf1SME3D3GOSrbjKBlXp1zAKXRZixFdGgwPdR7gsovlOSZYRvjp8z6HkMFazNcHzDLR7wxf2AnvQpJpUvssrPSz96xRF37PmvdIwIlebYRemHwAaNRT_4zC-CFczGncHctZ4WThw0TrUNHpDT-3Z_IXzO1MG5y0wUBDxiF9GDoVoBRPthwXgYTgXugtHjt4Hw3DdeR7p-ATTmjbnGihYy0TPMJVxztb5mjx5_cxKDpS1aTllEtIQbDHfaayxl1nq3VPTts82w_TkhHvYyknqewe5c1eVZiUxYrNhJSnL1MfvCjaJOY3mLcYkAxuMSzkMY7zm3V73jEtIUlItGap3rTI_ZbC45vMxPQnF_CpI6EyBTcuxrvINgckLvYknpcFAzl2k6J3N_Mm5DjzswE080byrC3FSD-G_i
# load the log image rationing
# stack_after_before_log<-
# stack(here::here("MAXAR_data","image_differencing","stack_after_before_log.tif"))
#
# # Extract all values from the raster into a data frame
# rstDF <- values(stack_after_before_log)
#
# head(rstDF)
# # Check NA's in the data
# idx <- complete.cases(rstDF)
# head(idx)
#
# # Initiate the raster datasets that will hold all clustering solutions
# # from 2 groups/clusters up to 12
# rstKM <- raster(stack_after_before_log[[1]])
# rstCLARA <- raster(stack_after_before_log[[1]])
#
# # increase the memory limit to handle the processing
# memory.size(40000)
#
# for(nClust in 2:12){
#
# cat("-> Clustering data for nClust =",nClust,"......")
#
# # Perform K-means clustering
# km <- kmeans(rstDF[idx,], centers = nClust, iter.max = 50)
#
# # Perform CLARA's clustering (using manhattan distance)
# cla <- cluster::clara(rstDF[idx, ], k = nClust, metric = "manhattan")
#
# # Create a temporary integer vector for holding cluster numbers
# kmClust <- vector(mode = "integer", length = ncell(stack_after_before_log))
# claClust <- vector(mode = "integer", length = ncell(stack_after_before_log))
#
# # Generate the temporary clustering vector for K-means (keeps track of NA's)
# kmClust[!idx] <- NA
# kmClust[idx] <- km$cluster
#
# # Generate the temporary clustering vector for CLARA (keeps track of NA's too ;-)
# claClust[!idx] <- NA
# claClust[idx] <- cla$clustering
#
# # Create a temporary raster for holding the new clustering solution
# # K-means
# tmpRstKM <- raster(stack_after_before_log[[1]])
# # CLARA
# tmpRstCLARA <- raster(stack_after_before_log[[1]])
#
# # Set raster values with the cluster vector
# # K-means
# values(tmpRstKM) <- kmClust
# # CLARA
# values(tmpRstCLARA) <- claClust
#
# # Stack the temporary rasters onto the final ones
# if(nClust==2){
# rstKM <- tmpRstKM
# rstCLARA <- tmpRstCLARA
# }else{
# rstKM <- stack(rstKM, tmpRstKM)
# rstCLARA <- stack(rstCLARA, tmpRstCLARA)
# }
#
# cat(" done!\n\n")
# }
#
#
# # Write the clustering solutions for each algorithm
#
# save_kmean<-
# writeRaster(rstKM,
# filename=file.path(here::here("MAXAR_data","image_differencing"),"raster_kmean.tif"))
# save_clara<-
# writeRaster(rstCLARA,
# filename=file.path(here::here("MAXAR_data","image_differencing"),"raster_clara.tif"))
#
#
# #Evaluating Unsupervised Classification/Clustering Performance
#
# library(clusterCrit)
#
# # Start a data frame that will store all silhouette values
# # for k-means and CLARA
#
# clustPerfSI <- data.frame(nClust = 2:12, SI_KM = NA, SI_CLARA = NA)
#
# for(i in 1:nlayers(rstKM)){ # Iterate through each layer
#
# cat("-> Evaluating clustering performance for nClust =",(2:12)[i],"......")
#
# # Extract random cell samples stratified by cluster
# cellIdx_RstKM <- sampleStratified(rstKM[[i]], size = 2000)
# cellIdx_rstCLARA <- sampleStratified(rstCLARA[[i]], size = 2000)
#
# # Get cell values from the Stratified Random Sample from the raster
# # data frame object (rstDF)
# rstDFStRS_KM <- rstDF[cellIdx_RstKM[,1], ]
# rstDFStRS_CLARA <- rstDF[cellIdx_rstCLARA[,1], ]
#
# # Make sure all columns are numeric (intCriteria function is picky on this)
# rstDFStRS_KM[] <- sapply(rstDFStRS_KM, as.numeric)
# rstDFStRS_CLARA[] <- sapply(rstDFStRS_CLARA, as.numeric)
#
# # Compute the sample-based Silhouette index for:
#
# # K-means
# clCritKM <- intCriteria(traj = rstDFStRS_KM,
# part = as.integer(cellIdx_RstKM[,2]),
# crit = "Silhouette")
# # and CLARA
# clCritCLARA <- intCriteria(traj = rstDFStRS_CLARA,
# part = as.integer(cellIdx_rstCLARA[,2]),
# crit = "Silhouette")
#
# # Write the silhouette index value to clustPerfSI data frame holding
# # all results
# clustPerfSI[i, "SI_KM"] <- clCritKM[[1]][1]
# clustPerfSI[i, "SI_CLARA"] <- clCritCLARA[[1]][1]
#
# cat(" done!\n\n")
#
# }
#
# # save the results in a csv file
#
# save_silhouette<-
# write.csv(clustPerfSI, file = here::here("MAXAR_data","image_differencing","silhouette_clustPerfSI.csv"), row.names = FALSE)
#
# # print out a table with the silhouette index results for comparing each clustering solution:
#
# silhouette_table<-
# knitr::kable(clustPerfSI,
# digits = 3, align = "c",
# col.names = c("#clusters","Avg. Silhouette (k-means)","Avg. Silhouette (CLARA)"))
#
# silhouette_table
#
# # make a plot for comparing the two algorithms
#
# plot(clustPerfSI[,1], clustPerfSI[,2],
# xlim = c(1,13), ylim = range(clustPerfSI[,2:3]), type = "n",
# ylab="Avg. Silhouette Index", xlab="# of clusters",
# main="Silhouette index by # of clusters")
#
# # Plot Avg Silhouette values across # of clusters for K-means
# lines(clustPerfSI[,1], clustPerfSI[,2], col="red")
# # Plot Avg Silhouette values across # of clusters for CLARA
# lines(clustPerfSI[,1], clustPerfSI[,3], col="blue")
#
# # Grid lines
# abline(v = 1:13, lty=2, col="light grey")
# abline(h = seq(0.30,0.44,0.02), lty=2, col="light grey")
#
# legend("topright", legend=c("K-means","CLARA"), col=c("red","blue"), lty=1, lwd=1)
I won’t run the code here again as it takes a long time to finish, outcome again is saved.
Below is the silhouette analysis where we evaluate the performance of both algorithms
# get the silhouette coerfficient saved in the table earlier
clustPerfSI<-
read.csv(here::here("MAXAR_data","image_differencing","silhouette_clustPerfSI.csv"))
#clustPerfSI
#call the table where the solhouette coefficeitn are stored
knitr::kable(clustPerfSI, digits = 3, align = "c",
col.names = c("#clusters","Avg. Silhouette (k-means)","Avg. Silhouette (CLARA)"))
| #clusters | Avg. Silhouette (k-means) | Avg. Silhouette (CLARA) |
|---|---|---|
| 2 | 0.470 | 0.460 |
| 3 | 0.490 | 0.438 |
| 4 | 0.424 | 0.398 |
| 5 | 0.423 | 0.368 |
| 6 | 0.391 | 0.360 |
| 7 | 0.365 | 0.355 |
| 8 | 0.355 | 0.342 |
| 9 | 0.281 | 0.293 |
| 10 | 0.307 | 0.295 |
| 11 | 0.312 | 0.285 |
| 12 | 0.270 | 0.266 |
# set the layout
plot(clustPerfSI[,1], clustPerfSI[,2],
xlim = c(1,13), ylim = range(clustPerfSI[,2:3]), type = "n",
ylab="Avg. Silhouette Coeff.", xlab="Number of clusters",
main="Silhouette Coeff. by number of clusters")
# Plot Avg Silhouette values across # of clusters for K-means
lines(clustPerfSI[,1], clustPerfSI[,2], col="navyblue",lwd=2)
# Plot Avg Silhouette values across # of clusters for CLARA
lines(clustPerfSI[,1], clustPerfSI[,3], col="lightsteelblue4",lwd=2)
# Grid lines
abline(v = 1:13, lty=2, col="light grey")
abline(h = seq(0.30,0.44,0.02), lty=2, col="light grey")
legend("topright", legend=c("K-means","CLARA"), col=c("navyblue","lightsteelblue4"), lty=1, lwd=2)
We can observe that Kmean cluster 3 had higher silhouette coefficient, therefore, 3 classifications will be elaborated.
let’s plot the clustered result with 3 clusters.
the classification was interpreted through projecting the 3 clusters stak in QGIS on a real photo as per below. for example, severe damage can be easily classified as the port area damages were visible.
knitr::include_graphics(here::here("MAXAR_data","image_differencing","kmeans_classification vs. real image.png"))